#	  Author: Arturas Rozenas, ar71@duke.edu
#	  Date: Aug 2011
#	  Replication code for paper "Statistical Model for Party 	  System's Analysis", Political Analysis
#    THIS CODE REPLICATES (up to mcmc error) TABLE 2 and     FIGURES 2 and 3 in the paper
#    setwd to the current file folder

source("orcomp.R")
source("aux_functions.R")
load("cses_data") ## the data are from CSES (all three 						modules)
					# this is .R data list
					# contains: 
					# v - vote-share matrix
					# data - data variables
					
# CREATE DESIGN MATRIXES

v=data$v
X=cbind(1,log(data$m))
colnames(X)=c("I","log(m)")
n=apply(data$v,1,function(x)(sum(x>0)))
G=cbind(1,log(data$m),n)
colnames(G)=c(colnames(X),"n")

l.bound=1 ### there should be at least 1 party in a country (a safe assumption that makes sure we don't predict zero parties)

library(mnormt)
library(TeachingDemos)

n.sample=2000
ptm <- proc.time()
out=orcomp(n,X,v,G,l.bound=l.bound,n.sample=n.sample,burn=5000,thin=5)
proc.time() - ptm

#========================================================
# run the following to reproduce table 3
#========================================================

D=L=max(n)-1; K=ncol(X); N=nrow(X); l.bound=1

bpost=apply(out$b,2,mean)
gpost=apply(out$g,c(1,2),mean)

cis=function(x){emp.hpd(x,conf=.95)}

conf=t(apply(out$b,2,cis))
R1=cbind(bpost,paste("[",conf[,1],", ",conf[,2],"]",sep=""))
rownames(R1)=c(colnames(X),"ln(w)")

conf=t(apply(out$mu,2,cis))

R2=cbind(apply(out$mu,2,mean),paste("[",conf[,1],", ",conf[,2],"]",sep=""))
R=cbind(R1,R2)

v.pred=function(b,g,X){ ### predicted vote-share function
n.pred=n.med(b,X)
G.pred=c(X,n.pred)
y.pred=(G.pred%*%g)[1:(n.pred-1)]
return(c(y.v(y.pred),rep(0,L+1-n.pred)))
}

n.med=function(b,X){
min(round(median(rtnegbin(10^4,mu=exp(X%*%b[1:K]),dispersion=exp(b[K+1]),l.bound = l.bound))),D+1)
} #### evaluate median predicted count; cannot exceed the dimension of the model


dist=function(x,y){ ##### distance function
return(.5*sum(abs(x-y)))
}


V=matrix(NA,N,L+1) ###### in sample predicted voteshares
for (i in 1:nrow(X)){
V[i,]=v.pred(bpost,gpost,X[i,])
}


print(R)
cat("nu",mean(out$nu))
cat("Correctly predicted", mean(abs(V-data$v)<.05))
cat("MAE", mean(mean(abs(V-data$v))))


#==============================
#	Reproduce Figure 4
#==============================



v.star=function(m){
t=sample(1:n.sample,1)
return(v.pred(out$b[t,],out$g[,,t],c(1,log(m))))
}


m.star=c(1,25,50,100)
e=.1
pd=matrix(NA,n.sample,3)
for (i in 1:nrow(pd)){
pd[i,1]=ifelse(dist(v.star(m.star[1]),v.star(m.star[2]))<e,0,1)
pd[i,2]=ifelse(dist(v.star(m.star[3]),v.star(m.star[4]))<e,0,1)
pd[i,3]=ifelse(dist(v.star(m.star[1]),v.star(m.star[4]))<e,0,1)
}

C=signif(round(apply(pd,2,mean,na.rm=TRUE),2),2)
C

par(mfrow=c(2,2))
par(mar=c(2.5,4.5,2,3))
par(oma=c(.01,.01,.01,.01))
for (i in 1:length(m.star)){
vpred=v.pred(bpost,gpost,c(1,log(m.star[i])))
names(vpred)=c(1:(L+1))
barplot(vpred[1:7], main=paste("M",m.star[i],sep="="),ylim=c(0,.4),cex.lab=1.5)
}

x=1:3


x[1]=expression(paste(pi[epsilon],' = ',"0.54"))
x[2]=expression(paste(pi[epsilon],' = ',"0.09"))
x[3]=expression(paste(pi[epsilon],' = ',"0.83"))


par(fig=c(0,1,0,1), new=TRUE)
plot(1, xlim=c(0,1),ylim=c(0,1), type="n", axes=FALSE, xlab="", ylab="")

text(.4,.9,x[1],cex=1.3)
arrows(.25,.85,.52,.85,lwd=2,code=2,col="grey70")
arrows(.25,.85,.52,.85,lwd=2,code=1,col="grey70")

text(.4,.17,x[2],cex=1.3)
arrows(.25,.13,.52,.13,lwd=2,code=2,col="grey70")
arrows(.25,.13,.52,.13,lwd=2,code=1,col="grey70")

text(.44,.42,x[3],cex=1.3)
arrows(.3,.5,.53,.25,lwd=2,code=2,col="grey70")
arrows(.3,.5,.53,.25,lwd=2,code=1,col="grey70")